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Abstract 

By means of Monte Carlo simulations of extensive air showers (EAS), we have 
performed a comprehensive study of the shower to shower fluctuations affecting the 
longitudinal and lateral development of EAS. We split the fluctuations into physical 
fluctuations and those induced by the thinning procedure customarily applied to 
simulate showers at EeV energies and above. We study the influence of thinning 
on the calculation of the shower to shower fluctuations in the simulations. For 
thinning levels larger than iZthin = 10~^ — 10~^, the determination of the shower 
to shower fluctuations is hampered by the artificial fluctuations induced by the 
thinning procedure. However, we show that shower to shower fluctuations can still 
be approximately estimated, and we provide expressions to calculate them. The 
influence of fluctuations of the depth of first interaction on the determination of 
shower to shower fluctuations is also addressed. 

Key words: Cosmic rays. Extensive air showers. Ground detector. Simulation, 
Muon component. Electromagnetic component 
PACS: 96.50.S, 96.50.sd, 13.85.Tp 



1 Introduction 



Extensive air showers (EAS) have been studied over the last 70 years [T]. They 
result from the interaction in the atmosphere of high-energy protons and nuclei 
arriving from space. The product of these collisions are a set of secondary 
particles carrying a fraction of the primary energy. These secondaries move 
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through the atmosphere and interact again generating new secondaries. The 
process continues, increasing the number of secondary particles, until their 
energies are too low to contribute to the generation of new particles. Particles 
reaching ground are sampled with arrays of detectors, and their properties are 
used to infer the properties of the primary initiating the shower. Measurements 
of the electron and muon density, of the arrival time of the particles at ground, 
and of the depth at which the shower has the maximum number of particles 
(Xmax); givc information on the arrival direction, primary energy, and on the 
mass of the primaries [T] . 

The complexity of the cascade phenomena, and the poor knowledge of the 
hadronic interactions at very high energy [2], make the experimental deter- 
mination of the properties of the primaries very difficult. Moreover, primary 
particles with the same energy, mass and direction produce secondary parti- 
cles with parameters that vary from shower to shower. This feature is called 
"shower to shower fluctuations". An understanding of the shower to shower 
fluctuations will help to improve the interpretation of cosmic-ray data. 

The calculation of shower to shower fluctuations can in principle be addressed 
with Monte Carlo simulations of extensive air showers. However, the number 
of particles that are produced in an air shower at ultra high energy (above 
~ 10^^ eV) is so large (~ 10^°), that it is almost impossible to follow the propa- 
gation to ground level of all the secondaries in the Monte Carlo in a reasonable 
amount of time, or even to store the large amount of information produced. 
For this reason, a statistical sampling procedure called "thinning" [3] is used 
in the simulations. Thinning algorithms typically consist on propagating only 
a small, representative fraction of the total number of particles in the shower, 
assigning statistical weights to the sampled particles to compensate for the re- 
jected ones. However, thinning algorithms introduce artificial fiuctuations in 
the simulated showers, hampering the determination of the intrinsic, physical 
shower to shower fiuctuations with Monte Carlo simulations. For this reason 
the study of fiuctuations using Monte Carlo simulations is quite difficult and 
uncertain. This is of utmost importance in cosmic-ray physics, since an incor- 
rect assumption on the shower to shower fiuctuations can lead to systematic 
errors on the determination of the parameters of the primary particles. 

In this work we address the problem of determining the true, physical shower 
to shower fiuctuations in Monte Carlo simulations, and quantify the effect of 
thinning on their determination. We give expressions that allow the estimation 
of physical fluctuations from Monte Carlo simulations, even in the case of 
relatively strongly thinned showers. 

The paper is organized as follows: In Section [2] we describe the simulations 
performed in this work, and the thinning algorithm adopted. In Section [3] we 
identify the different sources of fluctuations in shower simulations. In Section H] 
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we perform a comprehensive study of the fluctuations in the longitudinal and 
lateral shower development, and give expressions that allow to separate phys- 
ical shower to shower fluctuations from the artificial fluctuations induced by 
the thinning procedure. In Section O we quantify the influence on the shower 
to shower fluctuations of the fluctuations of the depth of first interaction of 
the primary initiating the shower. Finally, we summarize our conclusions in 
section |6l In the Appendix we give an explicit mathematical derivation of the 
expressions presented in Section HJ 



2 The simulations 

In this work we have used the air shower simulation program, AIRES [^||5] . 
along with the hadronic model QGSJETOl [6J to simulate proton and iron- 
induced showers with primary energy 10^^ eV. As explained above, due to the 
large number of particles that are created in the simulation, AIRES includes a 
statistical sampling algorithm, that consists on propagating a small, represen- 
tative fraction of the total number of particles, assigning a statistical weight 
w to the sampled particles to compensate for the rejected ones. The weight is 
adjusted in such a way that both the total energy and the average number of 
particles is guaranteed to be conserved. 

Before the simulation starts, the user indicates, as an input to AIRES, the 
relative thinning level -Rthin- The thinning energy -Ethin - the energy below 
which the thinning process starts - is defined as £'thin=-Rthin x Ep where, Ep 
is the primary energy. For ultra high energy cosmic ray shower simulations 
convenient values for the relative thinning are -Rthin = 10~^ — 10~^, but the 
actual choice depends on the purpose of the simulation. The thinning level 
affects both the simulation CPU time and the size of the output produced 
in the simulation, both typically behaving linearly with R^^^^. If we increase 
-Rthin by a factor of 10 the simulation speeds up by a similar factor, and the 
output is reduced accordingly, but the price to pay is an enhancement of the 
artifical fluctuations in the simulated showers as discussed below. 

We describe here the thinning algorithm implemented in the AIRES code [S], 
originally due to Hillas [3]. At the beginning of the simulation, the primary 
particle is assigned a weight w = 1. Then the primary is propagated and inter- 
acts in the atmosphere producing n secondary particles. Before incorporating 
any secondary particle in the simulation, the energy of the primary Ep which 
has generated that secondary is compared to Ethm- If Ep > -Ethin, then all the 
secondaries with energy greater or equal than Etun are kept, and their weight 
is equal to the weight of the primary particle. Secondaries with energy less 
than -Ethin are kept with a probability pi = Ei/Etun (Ei is the energy of i^^ 
secondary), and their weight is adjusted so that Wi = (l/pi) x w, with w be- 
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ing the weight of the mother particle producing that secondary. On the other 
hand if Ep < -Ethm, it nieans that the particle came from a previous thinning 
operation. Then, one and only one of all the produced secondaries - say the j*^ 

n 

- is kept, with probability pj = Ej/ J2 Ei. Again, the weight of this particle is 

i=l 

increased by a factor Wj = (l/pj) x w. 

To avoid confusion between particles and weights, we identify an entry with a 
particle explicitely followed in the simulation which is associated a weight w. 
Hence, an entry represents w particles. It is important to stress that once the 
thinning energy is reached, the number of entries Ne is no longer increased in 
the shower processes (only one secondary particle is followed in each interac- 
tion), while the number of particles does however increase, since the weight 
of each entry typically increases in the showering process. When evaluating 
a physical observable, each entry must be weighted with its corresponding 
statistical weight. 

In AIRES, the thinning algorithm is complemented with an "extended thin- 
ning algorithm" [S], designed to ensure that all the statistical weights are 
always smaller than a certain positive number (other algorithms based on this 
same idea are possible, see for instance [7]). To ensure this, an external pa- 
rameter called statistical weight factor W is available in the simulation. To 
further optimize the procedure of sampling, separated weight factors for elec- 
tromagnetic (W(EM)) and heavy particles (W^(HADRONIC)), are defined. 
The parameter Vr(HADRONIC) is specified indirectly by the ratio: 



^ w(m 

ly(HADRONIC) ^ ' 



The default value of this ratio in the simulation is AEH = 88. The default 
value of the weight factor ly (EM) is 12. In this paper, we will use these default 
values unless otherwise specified. 

We have simulated proton and iron-induced showers with primary energy Ep = 
10^^ eV, zenith angle 6 = 0°, 30°, 45° and 60° and relative thinning -Rthm = 
10^^, 10~®, 10~^ and 10"*^ (the two latter -Rthin for proton only). We have also 
simulated showers with relative thinning of -Rthm = 10^'' but with ly(EM) = 
0.1 instead of the default value. Finally, we have also simulated two sets of 
proton and iron- induced showers at ^ = 0°, with fixed depth of first interaction, 
starting at the corresponding mean interaction depth for protons and iron at 
10^9 eV. 
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3 Fluctuations in EAS 



In a real shower or in a simulation of an EAS, there are a number of dif- 
ferent fluctuations that can occur. Rather generally, we can make a simple 
classification as shown in Table [U 

"Physical fluctuations" are those due to physical processes in the shower. Here 
we split them into those due to the first interaction, and those occuring in the 
secondary interactions, as is customary, and because it has recently been sug- 
gested that "universal" shower properties may emerge when considering only 
the fluctuations in the first interaction point |8]. Physical fluctuations occur- 
ing in the first interaction are further divided into those affecting the depth 
of the first interaction, and those that arise from fluctuations of multiplicity 
or inelasticity also in the first interaction. 



Physical fluctuations 


- Depth of first interaction. 




- Multiplicity, inelasticity, etc, in 1^* interaction. 




- Secondary interactions. 


Experimental fluctuations 


- Detector response. 




- Sampling fluctuations. 


Artificial fluctuations 


- Thinning. 




- Un-thinning. 



Table 1 

Classification of the fiuctuations in a shower, arising from the physical processes 
in the shower and the measurement process, and those that appear only in shower 
simulations. 

In the case of real data, fluctuations are enlarged due to the detector response, 
and to the fact that the detector usually only samples a small fraction of the 
shower front. This "sampling fluctuation" is a statistically well known problem, 
and sampling fluctuations are rather well studied [9|10j. We will not consider 
them in this work. Also the detector response introduces an additional source 
of fluctuations, which are detector dependent, and will not be considered here. 

On the other hand, Monte Carlo simulated data is affected by artificial fluc- 
tuations due to the thinning and un-thinning (re-sampling) procedures. For 
the purposes of this work we do not need to consider the effect of fluctuations 
induced by the unthinning procedure [9]. 



5 



4 Fluctuations of the longitudinal and lateral shower development 



4-1 Fluctuations of the longitudinal profile 



74 = zi 3. (3) 



In Fig. [T] we show the average longitudinal profile of the number of elec- 
trons (left panels) and muons (right panels) N, obtained in simulations of 100 
proton-induced showers with Ep = 10^^ eV, for thinning levels -Rthin = 10~^ 
and 10~^ and 6 = 0°, 60°. Also shown are the relative shower to shower fluc- 
tuations a/N for electrons and muons. As it is well known [llj, the relative 
fluctuation has a minimum close to the depth of shower maximum. Also and 
as it is apparent from the figures, the dependence of the relative fluctuation 
on the thinning is small, at least for depths close to the depth of maximum. 

In Fig.[2]we show for the same showers in Fig.[Tl the skewness and the kurtosis 
of the distribution of the number of particles at different depths. It is worth 
recalling that the skewness of the distribution of a variable x is defined as 

= <(^. (2) 

where x (cr^) is the average (standard deviation) of x. The kurtosis is defined 
as 

{{x — xY) 
~i 

where the "—3" in the definition is a convention to make 74 = for a Gaussian 
distribution. Both the kurtosis and the skewness can be positive or negative. 
The skewness is a measure of the asymmetry of the distribution with respect 
to the mean value. A negative sign implies that the distribution is "deformed" 
towards values of x smaller than the mean. The contrary applies for a positive 
sign. The kurtosis is a measure of the length of the tails of the distribution. 
Positive values imply that the distribution has tails longer than those of a 
Gaussian, while negative values imply that the tails are shorter (for instance 
a flat distribution, a box, has kurtosis -1.2). A Gaussian distribution has 73 = 
74 = 0. 

Several remarks can be made from Fig. |2l Firstly, it is apparent that both, 
the skewness and the kurtosis of the distribution of the number of particles 
depend strongly on the thinning level, contrary to what happens to the mean 
N and to the relative fluctuations a/N . Close to the depth of shower maximum 
both the skewness and the kurtosis have local extrema very different from 
zero, implying that the distribution is strongly non-Gaussian. The skewness is 
negative and this implies that the distribution is asymmetric towards smaller 
values of than average. The positive values of the kurtosis imply that the 
distribution of A^ has tails longer than those of a Gaussian, at least close to 
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shower maximum. Remarkably, the log-Gaussian distribution, widely used to 
parameterize fluctuations in the number of electrons, has both 73 > and 
74 > 0, while the fluctuations predicted by Monte Carlo simulations near 
the maximum of the shower, have negative skewness. For muons and at large 
depths the skewness is close to zero, so that a Gaussian or a log-Gaussian 
distribution is a good approximation. For electrons, this is never the case. 



4.2 Fluctuations of the lateral profile at ground 

Of special importance for cosmic-ray physics performed with arrays of detec- 
tors is the study of fluctuations in the number of particles at ground. 

In Fig. [3] we show the relative fluctuations (cx/iV) of the total number of 
electrons (left panels) and muons (right panels) at ground (upper panels), and 
in a ring of width Ar at a distance r = 1000 m from the shower axis (lower 
panels). In both cases the fluctuations are shown as a function of the number 
of showers simulated. The ring was taken from rmin = 912 m to r^ax = 1092 
m, i.e. Ar = 180 m corresponding to a symmetric interval in the logarithm of 
r around r = 1000 m, chosen so that it compensates the decreasing density of 
particles with a larger area as r increases. As expected, fluctuations in the ring 
Ar are larger than the fluctuations in the whole ground. Also, the fluctuations 
in the ring have a stronger dependence with the thinning level used than those 
in the whole ground. This is easy to understand. An entry of weight w falling in 
the ring represents w particles, so that by losing or gaining just a single entry, 
one would lose or gain w particles and the fluctuations are enlarged. This effect 
is not so strong when accounting for all the particles falling anywhere on the 
ground. In Fig. 12] it can also be clearly seen that no reliable evaluation of the 
fluctuations can be done with less than about 20 simulated showers, especially 
in the case of fluctuations in the ring. It can also be seen that a thinning level 
of -Rthin = 10~^ or larger, introduces large artificial fluctuations, so that the 
shower to shower physical fluctuations cannot be evaluated reliably. This is 
however no obstacle to approximately estimate the physical shower to shower 
fluctuations as will be shown in the following. 



4-2.1 Physical shower to shower fluctuations at ground 

In the Appendix we prove that the distribution of the number of particles 
as obtained in thinned Monte Carlo simulations of extensive air showers has 
a mean N and a standard deviation a given by: 

N = N,w, (4) 
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and, 
where: 



(5) 



• Ne and s are respectively the mean and the standard deviation of the dis- 
tribution of the number of entries falhng in a given ring around shower 
axis in each shower, i.e., the distribution of the number of non-thinned (ex- 
phcitely sampled) particles in the simulation. 

• w and Q are respectively the mean and the standard deviation of the dis- 
tribution of weights assigned to the entries. 

Of course Eq. is exact since the thinning algorithm is designed to reproduce 
it. For Eq. ([5]), the proof only assumes that the probability for an entry to 
have a given weight w is independent of the probability of a shower to have a 
given number of entries Nf.. This is only approximate since the total number 
of entries and their weights are constrained by energy conservation. 

One can interpret Eq. as follows. If all entries (sampled particles) had the 
same weight equal to w in all the simulated showers, then the distribution of 
weights would not fluctuate from shower to shower and we would get = 0. 
In this limit we would obviously have a = ws. This special case of Eq. ([5]) 
was also found in [12]. In the particular case in which all weights are equal 
to w = 1, i.e. the shower is fully simulated and the thinning procedure is not 
applied, then clearly N = N^, a = s and the fluctuations would be obviously 
dominated by the true, physical shower to shower fluctuations. In the opposite 
limit, if we imagine that showers always have the same number of entries equal 
to Ne, i.e. there are no physical shower to shower fluctuations, then we would 
get s = 0, and the fluctuation in the number of particles would be solely due 
to the fluctuations of the weight of the entries, i.e., the fluctuations would 
be dominated by the thinning procedure and o"^ = iVef2^. Therefore, we can 
identify the flrst term of Eq. ([5]) with the artiflcial fluctuations introduced by 
the thinning procedure, and the second term with the true (physical) shower 
to shower fluctuations and deflne: 

= m', (6) 

and 

^l^ys = ^'S\ (7) 

Eq. ^ allows us to split the fluctuations into artiflcial and true ones and 
hence to estimate the effect of the thinning in a particular set of simulations, 
performed even with a relatively large value of the thinning level -Rthin ~ 
10~^, without the need to run new, more time-consuming and sometimes even 
impractical simulations with smaller -Rthin- 
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In the following, and by means of our Monte Carlo simulations, we numerically 
verify two key elements. Firstly, that Eq. ([5]) accounts for all the fluctuations 
(artificial and physical) appearing in simulations of EAS with thinning; and 
secondly that as -Rthin 0, (and the effect of the thinning procedure on the 
fluctuations is decreasingly important so that the fluctuations are increasingly 
dominated by the physical shower to shower fluctuations), the second term 
in Eq. ([5]) suffices to describe the fluctuations obtained in the Monte Carlo 
simulations. 

To see this in detail, we have calculated in Monte Carlo simulations, the av- 
erage weight w and the sigma of the distribution of weights Q, the average 
number of entries Ne and the corresponding sigma of its distribution s, as well 
as the average number of particles N and the sigma of its distribution a, both 
for electrons and muons, and compared them to what is predicted by Eq. 1^. 

Firstly in Fig. H] we show the average number N of electrons (left panel) and 
muons (right panel) versus the distance to the shower axis r for different 
thinning levels. As can be seen, the result is rather independent of the thin- 
ning level. This is not the case for the relative fluctuations a/N shown in 
Fig. El which depend strongly on the thinning level used in the simulations. 
For both electron and muons a/N has a minimum at the distance at which 
the number of particles is largest. Also, as expected, the fluctuations can be 
seen to converge to a common value at each distance as the thinning level 
decreases, because the effect of thinning is increasingly less important. The 
artificial fluctuations introduced by thinning in the case of electrons, do not 
contribute equally to the total fluctuation at all distances from the core as 
expected. It can be seen for instance that the relative fluctuation rises with r 
but the increase is smaller the smaller the thinning level. 

In Fig. [H] we plot the average weight w assigned in the process of thinning to 
electrons (left panel) and muons (right panel) as a function of the distance 
to the shower axis. In both cases the average weight is simply proportional 
to the thinning level w oc -Rthin, as expected. The average weight of electrons 
is typically 50 times larger than that of the muons. For electrons the average 
weight decreases at large distances to the core, because far from the core the 
electrons are mainly produced by muon decay, and muons carry a smaller 
weight. For the muons, there is a mild increase with r because the highest 
energy muons, which typically carry a small weight (i.e. they are less thinned 
than lower energy muons), are typically produced close to shower axis. 

In Fig. [7] we plot the relative fluctuations of the distribution of weights Q/w, 
for electrons (left panel) and muons (right panel) as a function of distance to 
the shower core. As expected the fluctuation of the weight decreases as -Rthin 
decreases and the showers are less thinned. Also, for values of -Rthin < 10~^ 
the relative fluctuation Q/tv is roughly independent of the thinning level, and 



9 



as a consequence we have that approximately Q oc -Rthin- 

In Fig. [HI we show the average number of entries Ne for the same simulations as 
in Fig. E] above. Clearly, we have that oc -R^-n as imposed by the constraint 
in Eq. that the average number of particles N has to be independent of 
-Rthin, together with the fact that w oc -Rthin- In Fig- El we show the relative 
fluctuation in the number of entries s/Ne- For small thinning levels (-Rthin < 
10~^), we find that the relative fluctuation is approximately independent of 
the thinning level, implying that s oc -R^J^j^- 

Finally, in Fig. [101 we compare the relative fluctuation of the number of parti- 
cles a/ N obtained directly in Monte Carlo simulations, with that predicted by 
Eq. ([5]), using the values of f2, Ne and s obtained in the same simulations. 
The comparison is shown for thinning levels -Rthin = 10^^, 10~^ and 10"''. 
The agreement between the a/N obtained in Monte Carlo simulations, and 
that predicted by Eq. ([S]) is at the level of < 20% for electrons and < 5% for 
muons, confirming that Eq. ^ accounts for all the fluctuations (artificial and 
physical) appearing in the simulations of EAS with thinning. 

From the scalings with -Rthin of the different magnitudes involved in Eq. ([5]) 
obtained before, it is straightforward to deduce that cTthin oc -Rthin, while Cphys 
should be approximately independent of -Rthin- This is seen in Fig. [TT] For 
thinning levels -Rthin ~ 10~^ and smaller, (Tphys is almost independent of -Rthin, 
while (Tthin depends strongly on -Rthin- In Fig- [HI it can also be seen that as 
the thinning level decreases the cxphys term increasingly dominates. This is of 
course expected, but it is remarkable that it was obtained from Monte Carlo 
simulations, and therefore it gives a strong support to our identification of 
(jphys with the true physical fluctuations. 



4- -2. 2 Dependence of fluctuations at ground on the number of particles 

Let us now consider the dependence of the fluctuations on the size of the 
ring around a distance to the shower core r, where particles are collected 
in the simulation. Let us consider how the density of particles is evaluated. 
If N{r, Ar) is the average number of particles at distance r in a small bin 
Ar (here we assume cylindrical symmetry around the shower axis, but the 
argument does not depend on this simplification), then the density of particles 
p(r) can be defined as 

/ , , N(r,Ar) , , 

p(r) = limA.-.o ; \ ^ (8) 
ZTrrAr 

In the limit Ar — )■ 0, pir) is finite, at least for r 7^ 0. However, the same is not 
true for the fluctuations, so that in general one can not define a "density of 
fluctuations" Pa{r). We can see this in a simple example. Assume that a{r, Ar) 
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is the standard deviation of the distribution of the number of particles in a 
bin of size Ar and at a distance r. One could try to define the density of 
fluctuations as 

p.(r) = lim^.^o (9) 
InrAr 

If the fluctuations in the number of particles at ground were purely Poissonian, 
we would have 

a(r, Ar) = ^N{r,Ar), (10) 

and therefore, 

jN{r, Ar) J27irArp{r) 
P.{r) = hmA._o 2vrrAr = 2vrrAr = ^ 

(11) 

i.e., we can not define a density of fluctuations for Poissonian shower to shower 
fluctuations in the number of particles. The dependence on bin-size has been 
identified with the fractal structure of showers [T3j . In the case of Poissonian 
fluctuations in the number of particles, the deduced behaviour is Pa{r) oc 
Ar" with a = —1/2, and we would be tempted to identify the coefficient 
a with a fractal exponent. However notice that for a Poissonian process no 
fractal structure is implied at all and it would be erroneous to call it a fractal 
exponent. 

On the other hand, in the case in which the fluctuations behave as: 

a(r,Ar) = /(r)Ar + 0(Ar2), (12) 

where /(r) is a function that does not depend on Ar, one can define a density 
of fiuctuations as can be shown trivially applying Eq. (Q. Remarkably, this is 
precisely the case of Furry's fiuctuations in which a ^ N = 27rrp(r)Ar 
and then p^r would be independent of Ar, i.e. a = 0. Recall that Furry statistics 
appears as a extremely simplified model of shower fiuctuations [H] , but it does 
take into account the branching structure of the shower (and therefore has an 
implicit fractal structure included). 

To our knowledge, the actual behaviour of the fiuctuations of showers and its 
dependence with the bin size is an open theoretical problem, with the theoret- 
ical prejudice ranging between two extremes: purely Poissonian fiuctuations 
a ~ V^, and stronger fiuctuations a ^ N. For instance, for the longitudinal 
development of showers one can show that in fact both types of behaviour 
occur [15], i.e., 

= aN^ + bN, (13) 
where a and b vary slowly with primary energy [15]. Near the maximum of 
the shower, the first term dominates and fiuctuations are not Poissonian. For 
the lateral distribution no such result exists but one would expect a similar 
conclusion. 
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In Fig. [T2I we show the relative fluctuation in the number of particles cx/iV, 
as a function of the bin size Ar, for several zenith angles 9 and distances r 
to the shower axis. Note that a/N is equal to p^- in the limit of Ar — )■ 0. A 
fit to a power law dependence (as suggested by the discussion above) gives 
a/N (X Ar" with a = —1/2 for both electrons and muons. As explained 
above this is suggestive of the fluctuations being Poissonian. However, if we 
split the fluctuations with the aid of Eq. ([5]) into thinning fluctuations, and 
physical fluctuations, it can be seen in Fig. [13] that Uphys/A^ is consistent 
with being flat with Ar, whereas athm/N behaves as a power law (a fit gives 
o"thin/A^ oc Ar^^/^). These results suggest that the artificial fluctuations are 
Poissonian, while physical fluctuations behave as Uphys oc N. 



5 Dependence of shower to shower fluctuations on composition and 
depth of first interaction 

In this Section we study the influence of the fluctuations in the depth of 
first interaction on the overall shower to shower fluctuations of the number 
of particles. In Figs. [H] and [15] we plot the relative fluctuations a/N in 10^^ 
eV proton and iron-induced showers respectively. In all panels we show the 
results of our regular simulations, together with the results of a special set of 
simulations performed by fixing the depth of first interaction of the primary 
particle (proton or iron) at the value of its mean interaction depth predicted 
by the QGSJET model (namely 44.9 g/cm^ for proton at 10^^ eV and 10.7 
g/cm^ for iron at the same energy). In all cases we use Eq. ([5|) to split the 
fluctuations into artificial and physical fluctuations, and we also show them 
in the figures. 

Firstly, it is interesting to see that the artificial fluctuations in the number 
of electrons or muons in iron showers are approximately equal to those in 
proton showers, while the physical fluctuations are smaller in iron than in 
proton-induced showers. The latter observation is a well-known effect which 
is attributed to the fact that showers initiated by a nuclei can be considered, 
in a first approximation, as a superposition of A (atomic mass) nucleons, each 
with an energy E/A with E the energy of the primary nucleus. 

It is rather remarkable that the relative fluctuations in the number of particles 
a/N m the two different sets of simulations (fixing or varying the depth of 
the first interaction point), are essentially the same. This conclusion applies to 
both the number of electrons and the number of muons. One could think that 
this is due to the fluctuations induced by thinning which mask the effect of the 
fluctuations of the depth of first interaction, however this does not seem to be 
the case, since as can be seen in Figs. [Tl]and[T5l neither the first term of Eq. ([5|) 
(the thinning fluctuations), nor the second term (the physical fluctuations) 
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change much when varying or fixing the depth of first interaction. We conclude 
that the relative shower to shower fluctuations on the number of particles at 
ground a/N are rather insensitive to the physical fluctuations of the depth of 
flrst interaction. This is related to the fact that the maximum of a shower at 
^ = 0°, where the fluctuations are minimum [IT], occurs near the ground. For 
other zenith angles a small difference appears in the physical fluctuations of 
the simulations performed with flxed and fluctuated flrst interaction point. 



6 Conclusions 

In this work we have performed a comprehensive study of shower to shower 
fluctuations by means of Monte Carlo simulations of extensive air showers. An 
understanding of the shower to shower fluctuations will help to improve the 
interpretation of cosmic-ray data. 

We have shown that the determination of the true, physical shower to shower 
fluctuations is hampered by the thinning procedure necessary to simulate in 
a practical manner air showers at EeV energies and above. However, we also 
show that the artiflcial fluctuations induced by thinning (o"thin) can be identi- 
fied and splitted from the physical fluctuations (uphys) with the aid of Eq. ([5]), 
which we have shown to account for all, true and artiflcial fluctuations ap- 
pearing in the simulations. Eq. (E]) reproduces the expectation that as the 
thinning level decreases -Rthm 0, and showers are less thinned, then the 
artiflcial fluctuations decrease, the physical ones become dominant, and they 
do not depend on i?thm- 

Our simulations also indicate that the physical shower to shower fluctuations 
of the number of particles at ground behave proportionally to the number 
of particles iV, while the artiflcial fluctuations are Poissonian, i.e., behave as 

Besides, we have shown that the size of the relative fluctuations due to the 
depth at which the flrst interaction initiating the shower occurs, is smaller or 
of the same order as the fluctuations occuring in the subsequent secondary 
interactions in the shower. 



7 Appendix 

In this Appendix, we calculate the probability distribution of the number of 
particles iV in a given bin of r distance to shower core, or of energy, making 
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some simplifying assumptions. We want to calculate the probability distribu- 
tion of particles, possibly in a given bin of r or of energy. We will assume that 
the probability for an entry (a non-thinned particle) to have a weight Wi is 
given by Pw{'Wi). In addition the number of entries, N^. in a given shower is a 
random variable with probability distribution Pe{Ne). Our main simplifying 
assumption is the following: we will assume that and Pg are independent 
of each other. This assumption is only approximate, because in a shower sim- 
ulated with thinning both the entries and the weight assigned to each particle 
are controlled by the branching of the shower and therefore they must be re- 
lated. However, as we will see a posteriori, the approximation is good enough 
for our purposes here, and it serves to clarify the role of thinning. 

Under this approximation we can write the probability P{N) of having N 
particles as, 

P{N) = P,{N,) \fdwi--- dwN^ P^{wi) ■ ■ ■ PUwnJ S{w, + --- + wn^-N) 

(14) 

where the (5— function expresses the constraint that the sum of weights is equal 
to the total number of particles. 

In what follows we evaluate this expression first by making further assumptions 
about the shape of Pe and P^, and afterwards in the general case using the 
characteristic function, related to the probability distribution. The definitions 
of cumulants and of the characteristic function can be found in any text book 
on statistics, for instance [T6l. 



We start with the integral 
I{N,, N)= fdwi--- dwN^ PUwi) ■ ■ ■ PUwnJ S{wi + --- + wn,-N), (15) 



and introduce the Fourier representation for the delta function. 
I{N,,N) = j dwi- --dwN^ PUwi) ■ --PUwnJ ^ fdk e*fe(-i+ -+-^.-iv). 



(16) 



Changing the order of integration gives 

I{N,,N) = ^jdk e-^^^l dwi PUwi)e'''"' 



dwM PJwnM 



2ti 



dke 



-ikN 



dwPu,{w)e 



ikw 



(17) 



To further continue with the evaluation of P{N) we need to make additional 
approximations. We assume that P^ is a Gaussian distribution with average 
w and rms Q 
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where A = 1/ \/2-itQ'^ is the probabihty normahzation. Its Fourier transforma- 
tion is given by 

J dw e'""^ P^w) = e'^"" e-'^'^'/l (19) 

Therefore, 

/(iV„iV) = — fdk e-^^(^~^«-) e-^^^'^'/2^ (20) 



21T 

For large values of N^, the sum in Eq. f[T^ can be approximated by an integral 
P{N) = ^ Pe(iVe) /(iVe, N) ^ f dN, P{N,) I{N,, N) . (21) 

If we assume that Pe{Ne) is also a Gaussian with average iVg and standard 
deviation s and inserting Eq. fl20l) in Eq. fl2T|) gives 

P(iV) = — [dk e-*'^^ / rfiV, g-(7v.-iv.)V{2s^) g^fciv.^ ^-k^N.nV2_ ^22) 

The integral over Ne can be done analytically 

P{N) = ^ fdk e-'^^ ex^[iN,wk - ]-{Ne^'' + s^w^)k'' + 0(A;3) ]. (23) 
2tt J 2 

where we neglect in the exponential powers of k larger than 2, after applying 
the saddle point approximation. We arrive at the final expression 

p(iV) = _^^-{N-Nfn2a-)^ (24) 



where 



V27rcr- 



N = NeW, 
J2 tTt" o2 I -.2 „2 



(x" = A^e^^ + w's\ (25) 



Notice that the above expressions have the correct asymptotic behaviour. If 
the particles have no weight, w — )■ 1 and i7 — t- 0, then the average number 
of particles is equal to the average number of entries (non-thinned particles) 
N = Ne and a = s. On the extreme case of a strongly thinned shower in 
which all particles are grouped together in a single entry [Ng = 1, s = 0) then 
N = w and a = f2, as expected. 

The above result is general and valid for any probability distribution for P^ 
and Pe- The only requirement is the "factorization" property given in Eq. (IT^ . 
From Eq. ( 1T71) . we introduce the characteristic function for the probability 
distribution P^, 

PUk) = I dw PUw) e""^, (26) 
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which, in general can be written as 

= exp[ ikai - ^^k'^ + ■■■]= e'3'^^\ (27) 

where the coefficients of the expansion of g{k) are related to the cumulants of 
the distribution of Pyj. For instance ai = w, a2 = etc. Then Eq. ([T^ reads 



P[N) = — fdk e-'^^ I dN, P,{N,) e'^=5(fc). (28) 

We define 

Pe(g) = / dNe Pe{N,) e^«^^ = exp[ thq - ^62?' + ■■■], (29) 

where as before 61 = and 1)2 = s^. Then we get 

P{N) =^ I dk e-^-^ PMk)) 
Ztt J 

=1. I dk e-''"^ exp[ thg{k) - )-h2g{kf + ■■■]. 

Where the function g{}z) = kw + i/2 k'^Q'^ + ■ • ■ . Then after some algebra 

P(N) = — [dk e-'^^ exp[ ikwN, - -k^s^w^ + N,n^) + ■■■] (30) 
2tt J 2 

The coefficients of the expansion around k = are again the cumulants of the 
distribution P{N), therefore we simply read the result given above in Eq. 



But, as a bonus, we obtain also all the other cumulants. For instance for the 
skewness we obtain 

73 = ^ = 4r X (mgiVe + 3ws^n' + w^M,), (31) 

where is the third central moment ( 1713 = {{w — w)^)) ) of the weight 
distribution and Mg is the third central moment of the distribution of the 
number of entries. In the same way, one can easily obtain other cumulants 
from the above expressions. 
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Figure 1. Upper panels: Longitudinal development of the average number of elec- 
trons (left) and muons (right) as a function of the slanted depth. Lower panels: 
Relative fluctuations a/N as a function of the slant depth for electrons (left) and 
muons (right). In all panels: 100 proton-induced showers with Ep = 10^^ eV and 
^ = 0° and 60°, were simulated with relative thinning iithin = 10"^, 10"'' and 10~^. 
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Figure 2. Upper panels: Skewness of the distribution of number of electrons (left) 
and muons (right) as a function of the slant depth. Lower panels: Same as the upper 
panels for the kurtosis. In all panels: 100 proton-induced showers of Ep = 10^^ eV 
with 9 = 0° and 60° were simulated with relative thinning i?th = 10~^, 10~^ and 
10-8. 
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Figure 3. Relative fluctuations (a/N) in the number of electrons (left panels) and 
muons (right panels) versus the number of showers in the simulations. In the upper 
panels we show the relative fluctuation of the total number of particles at ground. 
In the lower panels we show the fluctuations of the number of particles falling in a 
ring Ar at ground around a distance to the shower core r = 1000 m. In all panels: 
proton-induced showers of energy Ep = 10^^ eV and ^ = 0° were simulated with 
relative thinning as indicated in the insets. 
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Figure 4. Average number N of electrons (left panel) and muons (right panel) at 
ground as a function of the logarithm of the distance to the core. 100 proton-showers 
of energy 10^^ eV and = 0° were simulated with different thinning levels as 
indicated in the insets. 
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Figure 5. The relative fluctuation a/N in the number of electrons (left panel) and 
muons (right panel) at ground, as a function of the logarithm of the distance to the 
shower core, for the same sets of shower simulations as in Fig. [H 



Thinning: 


10^ 


Thinning: 




Thinning: 


10^ 




Figure 6. Average weight w of the distribution of weights of electrons (left panel) 
and muons (right panel) at ground, as a function of the logarithm of the distance 
to the shower core, for sets of 100 proton-induced showers at 10^^ eV and 6 = 0°, 
simulated with different thinning levels as indicated in the insets. 
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Figure 7. Relative fluctuations of the distribution of weights (Q/w) of electrons (left 
panel) and muons (right panel) at ground as a function of the logarithm of the 
distance to the core, for the same sets of shower simulations as in Fig. [6l 
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Figure 8. Average number of entries (non-tliinned particles) Ng for electrons (left 
panel) and muons (right panel) at ground, as a function of the logarithm of the 
distance to the core, for the same sets of shower simulations as in Fig. [6l 
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Figure 9. Relative fluctuation (s/Ne) of the distribution of the number of entries 
for electrons (left panel) and muons (right panel) at ground, as a function of the 
logarithm of the distance to the core, for the same sets of shower simulations as in 
Fig. El 
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Figure 10. Relative fluctuations a /N oi the distribution of number of electron (left 
panel) and muons (right panel) at ground, as a function of the logarithm of the 
distance to the core, for the same sets of shower simulations as in Fig. [HI The a/N 
obtained in the simulations (MC) is compared to that predicted by Eq. 



22 



o,^.^/W thining = 10' 
■ o,^.^/W thining = 10'^ 
a NJhining = 10'^ 
/N thin ing = 10'^ 
''^^Z N thining = 10'^ 
thining = 10'^ 




c,^.^/W thining = 10' 
■ o,h,n^W thining = 10'^ 

o,^.^/W thining = 10'^ 
' "^p/,^/^ ^'' ^ 

' "^p/,,,/^ ^'' ^ 

o /N thining = 10'^ 



Figure 11. Relative physical Uphys/-^ and thinning fluctuations o"thin/-^ (as predicted 
by Eq. ([5])) of the distribution of number of electron (left panel) and muons (right 
panel) at ground, as a function of the logarithm of the distance to the core, for the 
same sets of shower simulations as in Fig. [UJ 
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Figure 12. Relative fluctuations a/N of the distribution of the number of electrons 
(left panel) and muons (right panel) in a ring of width Ar, centered at different r at 
ground, as a function of the logarithm of the size of the bin Ar. 100 proton-induced 
showers of 10^^ eV energy were simulated at different 9 and flxed thinning level 

i?thin = 10-*^ 
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Figure 13. Relative fluctuations a/N of the distribution of the number of electrons 
(left panel) and muons (right panel) in a ring of width Ar, centered at different r 
at ground, as a function of the logarithm of size of the bin Ar. The upper lines cor- 
respond to simulations performed with relative thinning level -Rthin = 10~^ and the 
lower lines are for i2thin = 10"'^. The a/N obtained in the simulations is compared 
to that predicted by Eq. ([S]). The two terms in Eq. ((S]) corresponding to fluctuations 
induced by thinning (Jthin and physical fluctutations (Jphys are also shown, see insets. 
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Figure 14. Relative fluctuations a/N of the distribution of number of electron (left 
panel) and muons (right panel) at ground, as a function of the logarithm of the 
distance to the core, for 10^^ eV proton showers with ^ = 0°, and relative thinning 
level iJthin = 10~^. We show (squares) the result of fixing the first interaction depth 
at 44.9 g/cm^ (mean interaction depth of 10^^ eV proton-air collisions predicted by 
the QGSJETOl model), and also the case in which the depth of first interaction 
fluctuates (triangles). The a/N obtained in the simulations is compared to that 
predicted by Eq. ([5]) . The two terms in Eq. ([5]) corresponding to fluctuations induced 
by thinning cJthin and physical fluctutations Cphys ai^e also shown in all cases, see 
insets. 
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Figure 15. Same as Fig. [T3] for iron- induced showers. In the simulations with flxed 
first interaction depth, that depth was chosen at 10.7 g/cm'^ (corresponding to the 
mean interaction depth of 10^^ eV iron-air collisions predicted by the QGSJETOl 
model). Same symbols and line types as in Fig. [TH 
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